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ABSTRACT 

Lattice formulation of Finite Baryon Density QCD is problematic from com- 
puter simulation point of view; it is well known that for light quark masses 
the reconstructed partition function fails to be positive in a wide region of 
parameter space. For large bare quark masses, instead, it is possible to ob- 
tain more sensible results; problems are still present but restricted to a small 
region. We present evidence for a saturation transition independent from the 
gauge coupling (3 and for a transition line that, starting from the tempera- 
ture critical point at /x = 0, moves towards smaller (3 with increasing // as 
expected from simplified phenomenological arguments. 
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1 Introduction 



From the point of view of computer simulation, lattice approach to non per- 
turbative aspects of quantum field theory is a mature technique; apart from 
few exceptions, well consolidated schemes of simulation do exist, something 
like a recipes book, that allow studies, for example, of the most interesting 
features of QCD. The progress in the results is quite slow, in view of the 
large computing power needed for realistic calculation, but the field appears 
well founded. 

The mentioned "few exceptions", however, concern very interesting prob- 
lems, as well. The most paradigmatic of these dark zones is the study of 
thermodynamic of QCD in presence of non-zero baryonic density, shortly 
Finite Density QCD. The standard way to include the effects of baryonic 
matter on QCD vacuum leads to complex action in Euclidean formulation 
and this prevents the use of standard simulation algorithms, based on the 
idea of importance sampling, defined through a positive definite density of 
probability, e.g. the exponential of minus the Euclidean action. This problem 
can be rephrased stating the impossibility of defining a Boltzmann weight 
for each field configuration: only calculating the partition function we can 
define correctly the observables and obtain sensible results for quantities of 
physical interest. 

Calculations of partition functions are not infrequent in lattice simula- 
tions I], but their nature of extensive quantities raises the problem of the 
feasibility of this type of calculation with limited statistics, as forced from 
finite computing power. In the following section we will argue that, although 
reliable evaluation of the partition function of fermions coupled lattice gauge 
theories at zero baryon density is possible and successful 0, the extension 
of such technique for Finite Density QCD appears out of reaching for any 
reasonable statistics, at least in a range of theory parameters: for some val- 
ues of the chemical potential /i the phase of the fermionic determinant can 
be estimated only averaging over 0(e v ) configurations. 

In two recent papers [^, |J we have presented our results for Finite Den- 
sity QCD, mainly working with light quarks, in the infinitely strong coupling 
limit. Our approach consisted in trying to approximate the correct partition 
function using the modulus of the determinant, as suggested by the analy- 
sis of the SU(3) linear chain model ||, 0. The emerging picture is quite 
disappointing: we are essentially unable to reproduce the results of ref. [[?]], 
obtained with the MDP (Monomer-Dimer- Polymer) technique. Moreover we 
did show as the simulation scheme, known as Glasgow method |J, and pos- 
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sibly any other method based on the Grand Canonical Partition Function 
(GCPF) approach [[J, produce essentially the same results of our calcula- 
tions, when analyzed in a way to avoid perverse numerical effects due to 
rounding errors ||. 

In spite of these disappointing aspects, there are regions in the parameters 
space of the theory, in particular at large bare quark masses, in which one can 
hope to obtain reliable results, of some interest from a methodological point 
of view. In fact, at large quark masses (and any /3), the interval of chemical 
potential where the contribution of the phase can not be appreciated shrinks. 

The rest of the work, consequently, is devoted to the investigation of large 
quark mass limit, in some sense a favorite laboratory in which numerical 
techniques can be tested. Monitoring the expectation value of the phase 
of the Dirac determinant we can distinguish the regions in the parameter 
space where our evaluation of the partition function of finite density QCD 
is (in principle) exact from the ones where we miss a possible contribution 
to Z. A coherent picture seems to emerge from our data: a saturation 
transition exists at all couplings and merges, in the scaling region, to the 
true deconfining critical line that, with respect to fi = 0, moves towards 
smaller f3 with increasing fi. 

In the next section we will give arguments to explain why the contribution 
of the phase can not be measured and will present, in the strong coupling 
limit, a quantitative check of the Grand Canonical formulation results using 
data obtained with different techniques. The third section is devoted to the 
exposition of our approach to simulations of finite density QCD at finite 
coupling which exploits the main advantage of the MFA approach, i.e. the 
free mobility in the (/3, fx) plane. In the fourth section we present results 
for fermionic and gluonic observables, discussing the fate of the deconfining 
phase transition, expected, on phenomenological ground, when one increases 
the baryon density. The analisys is complemented with informations coming 
from the partially solved infinite bare mass limit |T0 |. In the last section a 



final discussion of the most important results is done. 

2 The partition function of Finite Density 
QCD 

The Finite Density QCD partition function can be written as 

Z = J[dU]e-^ {u) detA{U,m q ,fi) (1) 
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where, using the staggered formulation, the fermionic matrix A takes the 
standard form |TT] 



Ajj = m q 5ij + - Vv^Wv^Sj^+c* -Ul{i- v)bj,i-v\ 

^ ^=1,2,3 

The contribution of modulus | det A| of Dirac determinant and its phase 0a 
can be separated as 

Z = Z||(e^ A )|| (2) 

where 

Z\\ = J [dU]e-^ u) \ det A{U,m q ,v)\ (3) 

is the partition function of the model with the modulus of the determinant 
(modulus QCD in the following), and 





det A| 


Jc/> 


J[dU]e-0 s o(u) 


|detA 


1 



(4) 



It is clear from eq. (|2|) that, in the thermodynamical limit, the theory 
defined by means of Z\\ is physically different from the original theory only 
when the expectation value of the cosine of the phase of fermion determi- 
nant is vanishing exponentially with the system volume. In the regions of 
parameter space where the aforementioned expectation value is not 0(e~ v ) 
modulus QCD is an equivalent formulation of Finite Density QCD i.e. indis- 
tinguishable in the thermodynamical limit. In the rest of parameter space, 
modulus QCD clearly overestimates the true QCD partition function. 

Let us try to better illustrate this concept looking at figure 1. It refers 
to infinitely strong coupling limit (3 = and V = 6 3 x 4. At fixed quark 
mass m q the partition function of the system is only dependent on chemical 
potential fi. If we plot the free energy versus /i we can extract the phase 
structure from the appearance of a singularity in (some derivative of) the 
curve. 

Two extreme limits are well known. At /i = we get the logarithm 
of the usual fermion determinant averaged over gauge field configurations 
with a flat distribution: an average of a well defined (real and positive) 
quantity that can be computed. On the other hand in the large fi limit only 
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the last term of Grand Canonical Partition Function (see later) survives: 
det A — > (1/2) 3V e 3VtJl and the free energy is a straight line with slope 3V. In 
this limit the (baryon) number density, defined as: 

is equal to 1, and we can say that we are in a saturation regime, with Pauli 
exclusion principle preventing further increase of baryon density. In these 
two limits modulus QCD is coincident with the true theory and deviations 
are possible only in the intermediate region. 

Starting from \i — 0, we can use the data of fig. 5 in ref. 0, regarding 
number density at m q = 0.1, in order to reconstruct the free energy of the 
true theory as seen from the MDP approach. This is shown in figure 1 as 
the dotted line. If we superimpose the results of modulus QCD (continuous 
line) we can easily identify three regions: 

• n < fix — 0.3, which defines the onset in modulus QCD, where the 
number density is essentially zero; 

• \i > \ii — 1.0 , the saturated region; 

• Hi < n < /i 2 , the region where modulus QCD grossly overestimates the 
free energy of true theory. 

(as stated in [Q, using Glasgow prescription for dealing with the complex 
determinant, we obtain, for the free energy, exactly the same results as in 
modulus QCD). 

In figure 2 we report, for the same lattice and quark mass, the difference 
between the free energy of modulus QCD and the estimation based on data 
of ref. 0]. Superimposed to that we plot the expectation value of {e l ^^)\\ 
at the same value of the parameters. It is evident that the intermediate 
region is where the phase term is vanishing within statistical errors. If we 
concentrate on a value of fi inside this region, for example \l — 0.7, and we 
plot the distributions of the phase and the (logarithm of the) modulus of 
fermion determinant of single field configurations, we can see (figure 3) that 
modulus distribution is behaved as expected, while the phase distribution is 
almost flat. 

These distributions have been computed using N ~ 2500 configurations 
of a 6 3 • 4 lattice. With this statistics we can hope to measure accurately the 
phase term (e^ A )ii only down to 0(l/\/N) (~ 0.02 for our runs), far from the 
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0(e~ v ) order needed in principle. Even with a statistics of some thousands 
of configurations, we can say nothing on free energy of true theory in the 
range \i\ < \i < that covers the region where the number density varies 
rapidly. This does not imply necessarily that the phase is relevant in this 
region: for example it could go to zero as e~ Vs , with Vs the spatial volume, 
being in this case at the same time irrelevant and non measurable! 

The situation becomes somewhat better if we move to large quark mass: 
the range ([ii, (/,%), where finite statistics effects prevent to obtain a sensible 
evaluation of free energy, becomes narrower (see later), thus allowing the 
study of the model in a wider parameters range. The same scenario holds at 
finite coupling too, allowing us to investigate a great part of the parameter 
space. 

3 Simulation Scheme 

In this section we will present the simulation scheme that we have used in 
our work. Our simulations are based on the GCPF (Gran Canonical Parti- 
tion Function) formalism with an MFA (Microcanonical Fermionic Average) 
inspired approach for intermediate coupling analysis. 

The GCPF formalism allows one to write the fermionic determinant as a 
polynomial in the fugacity z = e M : 

detA(C/;m,/x) = det(G + e^T + e^T*) = z 3V det{P{U; m) - z' 1 ) 

3-Lf 
n=-3£3 

where the propagator matrix P |J is 

P(U;m)= ( "_ G T T T Q \ 

in which G contains the spatial links and the mass term, T contains the 
forward temporal links and V = L s Lt is the lattice volume. 

Once fixed the quark mass m q , a complete diagonalization of P allows 
one to reconstruct, trough a recursion algorithm, the coefficients a n , hence 
det A, for any value of the chemical potential //. Due to the Z(L t ) symmetry 
of the eigenvalues of P it is possible to write P Lt in a block form and we only 
need to diagonalize a 6L^ x 6L^ matrix. 
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This general method has been implemented in the framework of an MFA 
1^1 inspired approach. The basic idea in MFA is the exploitation of the 
physical equivalence between the canonical and microcanonical formalism 
via the introduction of an explicit dependence on the pure gauge energy in 
the computation of the partition function. Indeed (HD can be written as: 

Z(P, fi,m)=J dEn(E)e- 6VI3E < S F eff {^ m q ) > E (6) 

where 

n{E) = J [dU]5{QVE - S g [U]) (7) 
is the density of states at fixed pure gauge energy E, and 



< b eff {fi,m q ) > E = (M 



is the average over gauge field configurations at fixed energy E of a suitable 
definition of effective fermionic action. 

For the calculation of < S^tf >e we proceed as follows: first, we choose a 
set of energies selected to cover the range of f3 we are interested in. Secondly, 
for all the energies in the set, we generate gauge field configurations using 
a pseudo-microcanonical code; the generation of gauge fields at fixed energy 
is not the costly part of the whole procedure, so we can well decorrelate 
the configurations used for measuring the Dirac operator. Then, a standard 
NAG routine is used in order to obtain the complete set of eigenvalues of 
the propagator matrix P. At this point we can reconstruct the fugacity 
expansion coefficients a n or, without any substantial additional computer 
cost, use the eigenvalues to explore the possibilities offered by alternative 
prescriptions for the fermionic effective action, i.e. evaluate the modulus of 
the determinant and hence Zn. At the end, we have the fermionic effective 
action evaluated at discrete energy values: a polynomial interpolation allows 
the reconstruction at arbitrary values of the energy E, in order to perform 
the numerical one-dimensional integration in (||) and obtain the partition 
function Z\\ fi, m) . 

In a previous work Q we have found evidence for numerical instabilities 
in the standard evaluation of coefficients a n , whose origin lies on the ordering 
of the eigenvalues of P as calculated by a standard diagonalization routine. 
A random eigenvalue arrangement, before the calculation of the coefficients 
a n , is necessary in order to control rounding effects. In the present work we 
have always used this procedure to calculate the GCPF expansion coefficient. 
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To conclude this section let us, briefly, summarize the usefulness of MFA. 
This algorithm does not require a separate fermionic simulation for each value 
of (3, as the standard HMC (Hybrid Monte Carlo) algorithms, thus allowing 
us to extend the analysis to the whole relevant (3 range without an additional 
computer cost. 

Moreover, the basic idea of MFA is to consider the fermionic determinant 
(or its absolute value) as an observable. So det A is not in the integration 
measure and one avoids, in principle, the problem of dealing with a complex 
quantity in the generation of configurations. We have seen, however, that 
an unaffordable (eventually exponential) statistics is necessary to calculate 
Z in some \x range so that this advantage on direct simulation schemes can 
not, in general, be exploited. 

4 Large quark mass results 

In this section we will present results in the large bare quark mass limit both 
in the strong and intermediate coupling QCD at finite density. We have 
performed simulations in a 4 3 x 4 lattice (10 masses m q = 1.0 — > 5.0) in the 
range of the chemical potential \i G [0.0,4.0] and, for intermediate coupling, 
of (3 e [4.0,6.0]. Some (low statistic) data for a 6 3 x 4 lattice will also be 
shown. 

Firstly we have located the range (/ii,/z 2 ) (as a function of (3 and m q ), 
in which the average (^) is, with available statistics, indistinguishable from 
zero (see fig. 2). In figure 4 we report the gap A/i = ^2 — ^1 versus the quark 
mass computed fixing the gauge coupling value f3 = 5.5, a value near to the 
temperature induced transition. From figure 4 it is evident the tendency of 
A/z to decrease as the quark mass m q is increased. 

We have calculated the partition function of QCD in the (3— \x plane. This 
calculation is in principle exact for \x ^ [/ii(/3), /i 2 (/3)] and is complemented 
with data of modulus QCD in the region where we are not able to measure 
the phase. All the data are presented in such a way to make evident when a 
possible contribution from the phase has been discarded. 

In figure 5-a we report dN{p)jd\i evaluated at m q = 1.5 and different 
values of the gauge coupling [3. We have chosen two values of (3: [3 = 5.3 for 
which the system is inside the confining phase, and (3 = 5.7 in the deconfined 
one. 

The figure shows a sharp peak for the derivative, moreover the position 
of the peak does not move with (3. The same quantity for the larger lattice is 
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reported in figure 5-b. We do not attempt to do a finite size scaling analisys 
but only note that the height of the peaks grows considerably suggesting 
a (saturation) transition at /i = /if and independent from (3. The same 
scenario holds at smaller m q with the only difference that the peaks becomes 
broader for masses smaller than ~ 1.0 signalling the well known phenomena 
of early onset in the number density , |§ . 

These findings are to be compared with the one obtained analytically at 
infinite m q for which the infinite bare mass QCD partition function factorizes 


Z([3, /i) = K([3, /i) • Z PG (P) ■ Z(P = 0, //) (9) 

with 1Z a irrelevant factor in the zero temperature thermodynamical limit 
(1Z — > 1 for V = L 4 — > oo) and Z(f3 = 0, fi) developing a first order saturation 
transition at T = 0. 

Due to the independence of /if from /3, we will investigate its dependence 
on the bare quark mass from the behavior of N(fi) at (3 = 0, where we 
are able to compare with other results (analytical as well as numerical). In 
figure 6 we report /xf(m g ): at m q < 2 we approach a linear dependence for 
/if in good agreement with the numerical results by MDP while, at larger 
masses, /if coincide asymptotically with the large mass limit log(2m 9 ) [jTOJ 
(as well as with the 1/3 of the nucleon mass |fT3|). 

Coming to the phase structure in the (3 — fi plane it is evident that the 
critical line at constant /i = /if can not be the only one. In order to sep- 
arate the confined phase from the deconfined one we need a transition line 
that, starting at the zero density finite temperature critical point, eventually 
merges with the saturation transition at smaller values of j3. This critical line 
is the relevant one from a physical point of view in the sense that here we 
can obtain finite quantities for physical observables when the lattice spacing 
goes to zero. In the infinite mass limit this line is vertical (at the critical 
coupling of the finite temperature pure gauge theory, see formula @) so we 
can expect that, for our large masses, it moves only slightly to smaller (3. If 
this is the case it is crucial the possibility to rely on data at fixed /i and any 
(3 to extract relevant informations. 

The saturation of number density is a pure lattice artifact (a saturated 
lattice corresponds, in the continuum limit, to a divergent number density). 
Therefore, in order to search for evidence of the transition expected on phe- 
nomenological grounds, it is necessary to restrict our analisys to values of fi 
where the discretization effects are smaller i.e. to a number density smaller 
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than the lattice half filling value 1/2. In this region of the parameter space 
the phase gives no contribution (see fig. 2) and we can rely on our results. 

We have studied the plaquette energy E{(3, /x), the Polyakov loop P{(3, /x) 
and the number density as a function of (3. 

In figure 7-a,7-b we report E(/3,fj) and dE(/3, [x)/df3, evaluated at bare 
quark mass m q = 1.8 and at different values of the chemical potential /x < /xf . 
In fig. 7-a we can clearly see a rapid variation of the observable for all the 
values of /x; for the /x — curve this happens in correspondence with the 
pseudo-temperature transition of zero density full QCD. The critical gauge 
coupling moves to smaller f3 as we increase /x. This phenomenon is also 
evident as a sharp peak in the figure of the derivative (fig 7-b). It is tempting 
to interpret this as an evidence of a temperature induced phase transition 
extending at non zero values of /x. Also the behavior of the Polyakov loop 
points in this direction: as can be seen in figure 8, P^{P) changes rapidly at 
values of (3 consistent with the ones obtained from the energy. The number 
density gives a less clear signal since it is forced to be a constant function of f3 
at /x = 0. Nevertheless we can see in figure 9 that the plot of this observable 
is still consistent with previous findings for the gluonic quantities. It is useful 
to remark that plotting the same quantities at fixed (3 as a function of /x we 
would be practically unable to see any signal. 

Signals for a developing discontinuity in all these observables rely on 
data in the region where the contribution of the phase is negligible but sim- 
ilar behavior is found at larger values of /x too (where we can only compute 
the observables of modulus QCD). All these findings are consistent with 
phenomenological expectations for the temperature-density QCD phase dia- 
gram where, increasing the baryonic density, the critical temperature of the 
deconfinement phase transition decreases. On the lattice this translates in 
a critical line that, starting at the zero density-finite temperature critical 
point, continues at smaller values of (3 for finite /x. 

To conclude our analysis we report in figure 10 the (/3, /x) phase diagram 
of the theory at m q = 1.8 and V = 4 4 . We can see two (critical) lines; the 
horizontal one is due to saturation, while the other should be the physical 
one. If the same scenario holds in the zero mass limit we can expect that, 
as the lattice spacing goes to zero, the latter critical line extends in all the 
scaling window eventually coinciding with the former in the zero temperature 
limit. 
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5 Conclusions 



In this paper we have studied Finite Density Lattice QCD by means of nu- 
merical simulations. As well known this approach, probably the only one 
able to tackle the non perturbative effects leading to quark-gluon plasma 
transition, suffers severe problems due to the lack of hermiticity of Dirac 
operator for a single realization of gauge fields. 

In the first part of the paper we have shown as, for small quark masses and 
strong coupling, any numerical algorithm based on the GCPF approach gives 
results different from what expected in the region where the contribution of 
the phase can not be evaluated. To our understanding only a statistic expo- 
nentially large with the system volume (and a consequently high accuracy in 
numerical calculations) can solve this problem. 

Moving to large quark mass region we meet a much better situation and 
a large part of the parameter space becomes accessible to numerical sim- 
ulations. We get, independently from the gauge coupling f3, a saturation 
transition at a chemical potential /if well compatible with the one predicted, 
in the strong coupling regime, by previous numerical and analytical analysis. 
The new result is the evidence of another transition line that connects the 
previous one to the second order critical point of the four flavor /j, = theory. 
This has to be regarded as the lattice counterpart of the transition line in 
the temperature-chemical potential plane that should separate the standard 
hadronic phase from the quark-gluon plasma phase. 

For the first time we have got some evidences that the behavior of finite 
density lattice QCD can be consistent with standard phenomenological ex- 
pectations. Larger lattices could clarify the nature of the transitions but the 
volumes attainable with reasonable computer resources make this program 
not effective. To extend these results to the small quark mass region is im- 
possible since the contribution of the phase to the partition function becomes 
not measurable practically in the whole parameter space. At the end we have 
to conclude that, until now, finite density lattice QCD, far from providing 
non perturbative quantitative insights in the behavior of quarks and gluons, 
can at most give us some qualitative indication. 

This work has been partly supported through a CICYT (Spain) - INFN 
(Italy) collaboration. A.G. was supported by a Istituto Nazionale di Fisica 
Nucleare fellowship at the University of Zaragoza. 
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Figure Captions 



Fig. 1: Free energy at /3 — and m q = 0.1 of true QCD in the MDP 
approach (dotted line) and of Modulus QCD (continuous line). 

Fig. 2: Normalized difference between the free energy of Modulus QCD 
and the free energy of true QCD in the MDP approach (continuous 
line) superimposed to the expectation value of the determinant phase 
at (3 — and m q = 0.1. 

Fig. 3: Distributions of the logarithm of the modulus (a) and of the 
phase (b) of the fermionic determinant in a 6 3 x 4 lattice at (3 = 0, 
m q = 0.1 and \i = 0.7 with iV = 2500 configurations. 

Fig. 4: Width of the /i region (A/x) in which the QCD partition function 
fails to be positive versus the quark bare mass m q in a 4 3 x 4 lattice at 
= 5.5. 

Fig. 5: Derivative of the number density respect to the chemical po- 
tential in a 4 3 x 4 (a) and 6 3 x 4 (b) lattice at m q = 1.5 for f3 = 5.3 
(dashed line), and /3 — 5.7 (continuous line). 

Fig. 6: Saturation critical chemical potential /if versus the quark bare 
mass m q in a 4 3 x 4 lattice at /3 — 0. Dotted line is the large mass 
limit, the dashed one is the result of [[/J. 

Fig. 7: Plaquette energy E(/3,jj,) (a) and its derivative dE({3, fij/dp 
(b) evaluated in a 4 3 x 4 lattice at m q = 1.8 for /j = 0.0 — > 1.5 (from 
the right to the left) in steps of 0.1. Dashed line are for /j > jjli. 

Fig. 8: Polyakov loop P(/3,fi) evaluated in a 4 3 x 4 lattice at m q = 1.8 
for n = 0.0 — > 1.5 (from the right to the left) in steps of 0.1. Dashed 
line are for /x > H\. 

Fig. 9: Number density evaluated in a 4 3 x 4 lattice at m q = 1.8 for 
H = 0.0 — > 1.5 (from the right to the left) in steps of 0.1. Dashed line 
are for /j > jjli. 

Fig. 10: Phase diagram for the 4 3 x 4 lattice at m q = 1.8 in the (/3,/i) 
plane, the dotted line is for /j > 
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Free Energy 
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